DISCONTINUOUS GALERKIN METHODS FOR THE HELMHOLTZ EQUATION 

WITH LARGE WAVE NUMBER 
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Abstract. This paper develops and analyzes some interior penalty discontinuous Galerkin methods using piecewise 
linear polynomials for the Helmholtz equation with the first order absorbing boundary condition in the two and three 
dimensions. It is proved that the proposed discontinuous Galerkin methods are stable (hence well-posed) without any 
mesh constraint. For each fixed wave number k, optimal order (with respect to h) error estimate in the broken i?^-norm 
and sub-optimal order estimate in the L^-norm are derived without any mesh constraint. The latter estimate improves to 
optimal order when the mesh size h is restricted to the preasymptotic regime (i.e., k^h > 1). Numerical experiments are 
also presented to gauge the theoretical result and to numerically examine the pollution effect (with respect to k) in the 
error bounds. The novelties of the proposed interior penalty discontinuous Galerkin methods include: first, the methods 
penalize not only the jumps of the function values across the element edges but also the jumps of the normal and tangential 
derivatives; second, the penalty parameters are taken as complex numbers of positive imaginary parts so essentially and 
practically no constraint is imposed on the penalty parameters. Since the Helmholtz problem is a non-Hermitian and 
indefinite linear problem, as expected, the crucial and the most difficult part of the whole analysis is to establish the 
stability estimates (i.e., a priori estimates) for the numerical solutions. To the end, the cruxes of our analysis are to 
establish and to make use of a local version of the Rellich identity (for the Laplacian) and to mimic the stability analysis 
for the PDF solutions given in [23l [24l |35| . 

Key words. Flclmholtz equation, time harmonic waves, absorbing boundary conditions, discontinuous Galerkin meth- 
ods, error estimates 

AMS subject classifications. 65N12, 65N15, 65N30, 78A40 

1. Introduction. Wave is ubiquitous, it arises in many branches of science, engineer- 
ing and industry (cf. [22l [39] and the references therein). It is significant to geoscience, 
petroleum engineering, telecommunication, and defense industry. Mathematically, wave 
propagation problems are described by hyperbolic partial differential equations, and the 
progress of solving wave-related application problems has largely depended on the progress 
of developing effective methods and algorithms to solve their governing partial differen- 
tial equations. Among many wave-related application problems, those dealing with high 
frequencies (or large wave numbers) wave propagation are most difficult to solve numeri- 
cally because they are strongly indefinite and non-Hermitian and their solutions are very 
oscillatory. These properties in turn make it very difficult to construct stable numerical 
schemes under practical mesh constraints. Furthermore, high frequency (or large wave 
number) requires to use very fine meshes in order to resolve highly oscillatory waves, and 
the use of fine meshes inevitably gives rise huge, strongly indefinite, ill-conditioned, and 
non-Hermitian (algebraic) systems to solve. All these difficulties associated with high 
frequency (large wave number) wave computation still remain to be resolved and are 
not mathematically well understood in high dimensions (cf. [49]), although considerable 
amount of progresses have been made in the past thirty years (cf. [29l [36l |39] and the 
references therein). 

The simplest prototype wave scattering problem is the following acoustic scattering 
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problem (with time dependence e'^*): 

(1.1) -Au-k\ = f mR'^\D, 

(1.2) u = on dD, 

(1.3) ^( dr ^ '^^"^ ~ """"V ° as r = |a;| ^ cx), 

where Z) C R'^ (rf = 2, 3), a bounded Lipschitz domain, denotes the scatterer and i = v^— 1 
denotes the imaginary unit. A; G R is a given positive (large) number and known as the 



wave number, m™'^ is the incident wave. Equation (1.1 ) is the well-known Helmholtz equa- 



tion and condition (1.3) is the Sommerfeld radiation condition at the infinity. Boundary 
condition (1.2) implies that the scatterer is sound-soft. 

To compute the solution of the above problem, due to (finite) memory and speed limi- 
tations of computers, one needs first to formulate the problem as a finite domain problem. 
Two major approaches have been developed for the task in the past thirty years. The first 
approach is boundary integral methods (cf. [38] and the references therein) and the other 
one is artificial boundary condition methods (cf. [2H1II3] and the references therein). In the 
case of boundary integral methods, one converts the original differential equation into a 
(complicate) boundary integral equation on the boundary dD of the scatterer D. Clearly, 
the trade-off is that the original simple differential equation could not be retained in the 
conversion. On the other hand, artificial boundary condition methods solve the given 
differential equation on a truncated computational domain by imposing suitable artifi- 
cial boundary conditions on the outer boundary of the computational domain. Artificial 
boundary condition methods can be divided into two groups. One group of the methods 
use sharp artificial boundaries (i.e., the boundary has zero width), appropriate boundary 
conditions, which are called absorbing boundary conditions" , then are imposed on the 
boundaries (cf. pSj lHO]). The second group of artificial boundary condition methods allow 
the artificial boundaries to have non-zero width, such fatten boundaries are called absorb- 
ing layers, where the artificial boundary conditions are usually constructed in the form of 
differential equations which replace the original wave equations in the absorbing layers. 
The methods of this second group are called "PML (perfectly matched layer) methods" 
(cf. fT3]). In this paper, we shall adopt the absorbing boundary condition approach for 



problem (1.1)-(1.3). Extension of the results of this paper to the PML formulations will 
be given elsewhere. 

Specifically, in this paper we consider the following Helmholtz problem: 

(1.4) -Au-k\ = f mn:=ni\D, 

(1.5) u = on Yd, 

(1.6) ■p^ + iku = g ouTr, 

onn 

where {D c) fii C R^ = 2,3 is a polygonal/polyhedral domain, which is often taken 
as a d-rectangle in applications. Tr := dfli, Td = dD, hence, dfl = TrUTd- uq denotes 
the unit outward normal to dQ. The Robin boundary condition (1.6) is known as the 
first order absorbing boundary condition (cf. [28j). We remark that the case D = also 
arises in applications either as a consequence of frequency domain treatment of waves or 
when time-harmonic solutions of the scalar wave equation are sought (cf. flUi ITT]). 
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For many years, the finite element method (and other type methods) has been widely 
used to discretize the Helmholtz equation ( |1.4[ ) with various types of boundary conditions 
(cf. [IlElElElinilTniEIlEalETlEniEall^ reference therein). It is 

well known that in every coordinate direction, one must put some minimal number of grid 
points in each wave length I = 2Ti/k in order to resolve the wave, that is, the mesh size h 
must satisfy the constraint hk ^ 1. In practice, 6 — 10 grid points are used in a wave length, 
which is often referred as the "rule of thumb". However, this "rule of thumb" was proved 
rigorously not long ago by Babuska et al [37] only in the one- dimensional case (called the 
preasymptotic error analysis). The main difficulty of the analysis is caused by the strong 
indefiniteness of the Helmholtz equation which in turn makes it hard to establish stability 
estimates for the finite element solution under the "rule of thumb" mesh constraint. In 
[37] . Babuska et al also showed that the if^-error bound for the finite element solution 
contains a pollution term that is related to the loss of stability with large wave numbers. 
Later, Babuska et al addressed the question whether it is possible to reduce the pollution 
effect in a series of papers (cf. Chapter 4 of [36] and the reference therein). It should 
be noted that under the stronger mesh condition that k'^h is sufficiently small, optimal 
(with respect to h) and quasi-optimal (with respect to k) error estimates for finite element 
approximations of the Helmholtz problem were established early by Aziz and Kellogg in 
[5] and Douglas, et al in [261 EI] using the so-called Schatz argument [S] (also see Chapter 
5 of [H]), and a similar result was also obtained in [33] using an operator perturbation 
argument. 

The work of [3 EEl [37] shows that in the 1-d case, due to the pollution effect, the finite 
element solution for the Helmholtz problem (1.4)-(1.5) deteriorates as the wave number 
becomes large if the practical mesh condition fc/i < 1 is used. The situation in the high 
dimensions is expected to be the same (at least not better) although, to the best of our 
knowledge, no such a rigorous analysis is known in the literature. The detailed analysis of 
[2ni ET] also shows that the pollution effect is inherent in the finite element method and is 
caused by the deterioration of stability of the Helmholtz operator as the wave number k 
becomes large. In order to minimize or eliminate (if possible) the pollution and to obtain 
more stable and more accurate approximate solutions for Helmholtz-type problems with 
large wave numbers, various nonstandard and generalized Galerkin methods have been 
proposed lately in the literature. These methods can be categorized into two groups. The 
first group of methods use nonstandard or stabilized discrete variational forms to approx- 
imate the Helmholtz operator so that the resulted discrete problems have better stability 
properties. Methods in this group include Galerkin-least-squares finite element methods 
[ISl El] , quasi-stabilized finite element methods [TU] , and discontinuous Galerkin methods 
[U [TT] ES] • The second group of methods abandon the use of piecewise polynomial trial 
and test functions and replace them by global polynomials or non-polynomial functions. 
Methods in this group include spectral methods [15], generalized Galerkin/finite element 
methods [lOl |8], partition of unity finite element methods [H], and meshless methods 
[8]. We also note that another very different and intensively studied approach for high 
frequency wave computation is geometrical optics, which studies asymptotic (nonlinear) 
approximations of the Helmholtz equation obtained when the frequency (or wave num- 
ber) tends infinity. We refer the reader to [29] and the references therein for some recent 
developments in geometrical optics and its variants. 

The goal of this paper is to develop some interior penalty discontinuous Galerkin 
(IPDG) methods for problem (1.4)-(1.5) in high dimensions. The focus of the paper is to 
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establish the rigorous stabihty and error analysis, in particular, the preasymptotic error 
analysis. For the ease of presentation and to better present ideas, we confine ourselves 
to only consider the case of linear element in this paper. Such a restriction is also due to 
the consideration that we shall present /ip-discontinuous Galerkin methods for problem 
(1.4)-(1.5) in a forth coming paper [32] which extends the work of this paper to high 
order elements. Compared with existing DG methods for the Helmholtz equation in the 
literature, the novelties of our interior penalty discontinuous Galerkin methods include 
the following: First, our mesh-dependent sesquilinear forms penalize not only the jumps 
of the function values across the element edges/faces but also the jumps of the normal and 
tangential derivatives. Recall that penalizing the jumps of the normal (and tangential) 
derivatives helps but is not essential for the success of IPDG methods in the case of 
coercive elliptic problems (cf. [3l|25]), however, it contributes critically to the stability of 
the IPDG methods of this paper. Second, a small but vitally important idea of this paper 
is to take the penalty parameters as complex numbers of positive imaginary parts. This 
idea also contributes critically to the stability of the IPDG methods of this paper. As 
a result, essentially and practically no constraint is imposed on the penalty parameters. 
Since the Helmholtz problem is a non-Hermitian and an indefinite linear problem, as 
expected, the crucial and the most difficult part of the whole analysis is to establish the 
stability estimates (i.e., a priori estimates) for the numerical solutions. To the end, the 
cruxes of our analysis are to establish and to make use of a local version of the Rellich 
identity (for the Laplacian) and to mimic the stability analysis for the PDE solutions 
given in [231 |2ll |35] . Suppose f^i is star-shaped with respect to a point xoj. The key idea 
is to use the special test function Vuh ■ {x — xqJ (defined element- wise), which is a valid 
candidate for any IPDG method. We remark that the same technique was successfully 
employed by Shen and Wang in |33] to carry out the stability and error analysis for the 
spectral Galerkin approximation of the Helmholtz problem. 

In the past fifteen years, DG methods have received a lot attentions and undergone 
intensive studies by many people. As is well known now, DG methods have several advan- 
tages over other types of numerical methods. For example, the trial and test spaces are 
very easy to construct, they can naturally handle inhomogeneous boundary conditions and 
curved boundaries; they also allow the use of highly nonuniform and unstructured meshes, 
and have built-in parallelism which permits coarse-grain parallelization. In addition, the 
fact that the mass matrices are block diagonal is an attractive feature in the context of 
time- dependent problems, especially if explicit time discretizations are used. We refer to 
[31 SI [ni [191 [201 |25l lall ill ill ST] and the references therein for a detailed account on 
DG methods for coercive elliptic and parabolic problems. In addition to the advantages 
listed above, the results of this paper also demonstrate the fiexibility and effectiveness of 
DG methods for strongly indefinite problems, which was not well understood before. 

The remainder of this paper is organized as follows. In Section [2| notation and some 
preliminaries are described and cited. In particular, the sharp (with respect to k) stability 
constant estimates of [21] for the solution of problem (1.4)-(1.5) in high dimensions were 
recalled. These estimates are critical for obtaining explicit dependence of the error bounds 
on the wave number k. In Section [3| the IPDG methods of this paper are formulated. 
Both symmetric and non-symmetric IPDG methods are constructed. However, since the 
Helmholtz equation and its solution are complex-valued, the non-symmetric terms in the 
IPDG sesquilinear form do not cancel each other when two arguments of the form are 
taken to be the same function. Instead, their difference is a pure imaginary quantity. 
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This is a main difference between non-symmetric IPDG for coercive elliptic problems 
and for indefinite Helmholtz type problems. As a result, the penalty parameters need 
to be chosen as complex numbers with positive imaginary parts to ensure the stability 
in both symmetric and non-symmetric IPDG methods. Section |4] devotes to stability 
analysis for the IPDG methods proposed in Section [3| It is proved that the proposed 
IPDG methods are stable (hence well-posed) without any mesh constraint. In Section [sj 
using the stability result of Section |4] we derive optimal order (with respect to h) error 
estimate in the broken if^-norm and sub-optimal order estimate in the L^-norm without 
any mesh constraint. The latter estimate improves to optimal order when the mesh size h 
is restricted to the preasymptotic regime (i.e., k'^h > 1). In particular, for appropriately 
chosen penalty parameters, it is shown that the error in the broken if^-norm is bounded 
by Cikh + C2k^^^h^^^ if kh < 1. Numerical tests in Section [g] suggest that the error in the 
broken if^-norm may have a better bound Cikh + C2k^h? and it is possible to tune the 
penalty parameters to significantly reduce the pollution error. We note that in the case 
k'^h is sufficiently small, optimal order (with respect to h) error estimate in the broken H^- 
norm can be derived by using the Schatz argument as done in |3Hl 13, |26l ET] . In Section 
[6} we present some numerical experiments to gauge our theoretical error estimates, to 
numerically examine the pollution effect in the error bounds, and to test the performance 
of the proposed IPDG methods. 

2. Notation and Preliminaries. The standard space, norm and inner product notation 
are adopted. Their definitions can be found in |lll[T8l[TT]. In particular, (-, ■)q and (■, ■)s 
for S C dQ denote the L^-inner product on complex-valued L'^iQ) and -Z^^(S) spaces, 
respectively. (-, ■) := (-, ■)n and (-,■) := (-, Let 

Hl^{VL) := {u e H\n) : M = on F^} . 

Throughout the paper, C is used to denote a generic positive constant which is indepen- 
dent of h and k. We also use the shorthand notation A < i? and B > A for the inequality 
A < CB and B > CA. y4^i?isa shorthand notation for the statement A < B and 
B<A. 

We now recall the definition of star-shaped domains. 

Definition 2.1. Q C R"^ is said to he a star-shaped domain with respect to xq E Q 
if there exists a nonnegative constant cq such that 

(2.1) {x-XQ)-nQ>CQ MxedQ. 

Q C R'^ is said to he strictly star-shaped if cq is positive. 

Throughout this paper, we assume that Qi is a strictly star-shaped domain. In prac- 
tice, Qi is often taken as a (i-rectangle, which trivially is a strictly star-shaped domain. 
We also assume the scatterer D is a star-shaped domain with respect to the same point 
x^i as Qi does. This then implies that xq^^ ^ D G Qi. More precisely, we assume that 
there exist constants > and cd >0 such that 

(2.2) {x - xnj ■ nn > cn^ Vs G Fr and {x - xnj ■ no > cd Vx G F^?. 
Here hq and ud are the unit outward normals to the boundaries of Q and D, respectively. 



Under these assumptions the following stability estimates for problem (1.4)-(1.5) were 
proved in [23l IS! [35] . 

Theorem 2.2. Suppose Qi C R'' is a strictly star-shaped domain and D G Qi is a 
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star-shaped domain. Then the solution u to problem (1.4) -(1.5) satisfies 

1 



(2.3) 



+ ||5'l|L2(r 



for j = 0, 1 ifu G if^/^+^(fi) for some e > 0. (2.3) also holds for j = 2 if u e H'^{Vt). 

3. Formulation of discontinuous Galerkin methods. To formulate our IPDG methods, 
we first need to introduce some notation. Let he a family of triangulations of the 
domain Q := Qi \ D parameterized by /i > 0. For any triangle/tetrahedron K & Th, we 
define Hk '■= diam(ii"). Similarly, for each edge/face e of G Th, define he ■= diam(e). 
We assume that the elements of Th satisfy the minimal angle condition. We define 

set of all interior edges/faces of 
set of all boundary edges/faces of on Tr, 
set of all boundary edges/faces of on Yd, 
£h U £h = set of all boundary edges/faces of 7/j, 
Sl U 8^ = set of all edges/faces of except those on Tji. 

We also define the jump [v] of v on an interior edge/face e = dK fl dK' as 



cR 

cRD 
cID 



HI. 

If e G set Hie 



v\k 
v\k' 



v\k', 
- v\k, 



if the global label of K is bigger, 
if the global label of K' is bigger. 



v\e. The following convention is adopted in this paper 
WIe: ^ 



iv\K + v\K') ife = dKr]dK'. 



If e E S^, set {f }|e = v\e. For every e = dKndK' G S^, let Ue be the unit outward normal 
to edge/face e of the element K if the global label of K is bigger and of the element K' 
if the other way around. For every e G Sj^^, let = the unit outward normal to dQ. 

Now we define the "energy" space E and the sesquilinear form ah{-,-) on E x E as 
follows: 



E =11 H\K), 



(3.1) 
where 

(3.2) 
(3.3) 

(3.4) 



ah{u,v) := bh{u,v) + i[Jo{u,v) + Ji{u,v) + Li{u,v)) \/u,v G E, 



bh{u,v) 
Jo{u,v) 

Ji{u,v) 



To 
h. 



du 
drie. 



a 



\u\ 



dv 
drie. 



du 



drie 



dv 



(3.5) L,{u,v):=Y,Y. 



du 



dv 
d7i 
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and a is an /i-independent real number. 7o,e,7i,e, and (3i^e are nonnegative numbers to be 
specified later. {r|}^~J denote an orthogonal coordinate frame on the edge/face e G Eh, 



and 



Vm ■ tI stands for the tangential derivative of u in the direction r| . 
Remark 3.1. (a) Clearly, ah{-, ■) is a consistent discretization for —A since 



-Am, v) = 

ah{u,v) for all u G Hl^{9) fl H'^{9) and v e E. 

(b) If we regard ah{-,-) as a bilinear form on the subspace of real valued functions 
in Hyj^{Q), then ah{-,-) is symmetric when a = 1 and is non-symmetric when a ^ 1. 
In particular, a = —1 would correspond to the non-symmetric IPDG method studied in 

for coercive elliptic problems. In this paper, for the ease of presentation, we only 
consider the case a = 1, nevertheless the main results of the paper can also be extended 
to the case a 1. 

(c) The terms in i(^jQ(u,v) + Ji{u,v) + Li(u,v)) are so-called penalty terms. 

(d) The penalty parameters in \(^Jo{u,v) + Ji{u,v) +Li{u,v)) are i7o,e,i7i,e and i/3i,e; 
respectively. So they are pure imaginary numbers with positive imaginary parts. It turns 
out that if any of them is replaced by a complex number with positive imaginary part, 
the ideas of the paper still apply. Here we set their real parts to be zero partly because 
the terms from real parts do not help much (and do not cause any problem either) in 
our theoretical analysis and partly for the ease of presentation. On the other hand, our 



numerical experiments in Section 6. 5 indicate that using penalty parameters with nonzero 



real parts helps to reduce the pollution effect in the error. 

(e) Penalizing the jumps of normal derivatives (i.e., the Ji term above) for second 
order PDEs was used early by Douglas and Dupont [25] in the context of finite element 
methods, by Baker 01]/ (with a different weighting, also see fSl^ ) for fourth order PDEs, 
and by Arnold in the context of IPDG methods for second order parabolic PDEs. 
Arnold /<?/ also proposed and analyzed IPDG methods which penalize higher order normal 
derivatives. Note that we do not introduce boundary terms for e G Sj^ in Ji to ensure the 
consistency ofah{-,-) with —A. 

On the other hand, the idea of penalizing the jumps of tangential derivatives (i.e., the 
Li term above) seems is new. Later we will show that without Li term and Ji term in 
ah{-, ■) the IPDG methods of this paper are still stable and convergent but under a stinger 
mesh constraint, see Section\^and\^ 

(f) In this paper we consider the scattering problem with time dependence e"^*, that 
is, the signs before i's in the Sommerfeld radiation condition (1.3) and its first order 



approximation (1.6) are positive. If we consider the scattering problem with time depen- 



dence e that is, the signs before i's in (1.3) and (1.6) are negative, then the penalty 



parameters should be complex numbers with negative imaginary parts. 
Next, we introduce the following semi- norms/norms on the space E: 



(3.6) 



\l,h 



E 

Ken 



Wv 



|2 



(3.7) 



\l,h ■~ 



dv 
On. 



L2(e) 



d-1 



dv 

d?e 



L2(e) 
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(3. 



\i,h 



\i,h 



and 



he 
70,e 



dv 
drie 



L2(e) 



It is easy to see that 
dD = %. 

Clearly, the sesquilinear form au 



l,h 



(3.9) 



(3.10) 



Rea/i(f , v) 



|2 
\l,h 



are norms on E if dD ^ but only semi-norms if 

■) with (7 = 1 satisfies: For any v E E 

dv 



2Re (I 



lmah{v,v) = Jo{v,v) + Ji{v,v) + Li{v,v). 



With the help of the sesquilinear form a/i(-, ■) we now introduce the following weak 
formulation for (|L41)-(|L51): Find uE En Hl{n) n H^^^n) such that 



(3.11) ah{u,v)-k^{u,v) + ik{u,v)ri, = {f,v) + {g,v)^ 



yvEEnH'r^{n)nHU^). 



The above formulation is consistent with the boundary value problem (1.4)-( L6)) b ecause 
ah{-, •) is consistent with —A. It is clear that, if m G H'^{Q) is the solution of (1.4)-(1.5), 
then (3.11) holds for all v E E. 

For any K E Th, let Pi{K) denote the set of all linear polynomials on K. We define 
our IPDG approximation space V'^ as 



n ^iw- 

K&Th 



Clearly, V'' C E C L'^{n). But <^ H\n). 

We are now ready to define our IPDG methods based on the weak formulation (3.11 ): 
Find Uh E such that 



(3.12) ah{uh, Vh) - k'^iuh, Vh) + ik{uh, Vh)rR = (/, Vh) + {g, Vh)^, 



E V\ 



In the next two sections, we shall study the stability and error analysis for the above 
IPDG methods. Especially, we are interested in knowing how the stability constants and 
error constants depend on the wave number k (and mesh size h, of course) and what are 
the "optimal" relationship between mesh size h and the wave number k. 

4. Stability estimates. Since the Helmholtz operator is not a coercive elliptic operator, 
the stability estimates given in Theorem 



2.2 



for the solution of problem (1.4)-(1.6) is far 
from trivial. We refer the reader to [23l [2^ 135] for a detailed exposition in this direction. 
This difficulty is certainly inherited by any numerical discretization of problem (1.4)-(1.6). 
In fact, the situation usually is worse in the discrete case because piecewise polynomials (or 
piecewise smooth functions) are rigid, they are not as flexible as the PDE trial functions. 
On the other hand, DC approximation functions are much more flexible than Lagrange 
finite element functions because they do not have any continuity constraint, instead, the 
continuity of the numerical solutions is enforced weakly through mesh-dependent bilinear 
or sesquilinear or nonlinear forms. 

The goal of this section is to derive stability estimates (or a priori estimates) for 
scheme (3.12). To the end, momentarily, we assume solution Uh to ( 3.12[ ) exists and will 
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revisit the existence and uniqueness issues later at the end of the section. We hke to note 
that because its strong indefiniteness, unhke in the case of coercive eUiptic and parabohc 
problems (cf. [Sj HI [HI |25l EH HEl HSl H?]), the well-posedness of scheme (3.12 ) is far from 
obvious under practical mesh constraints. 



To derive stability estimates for scheme (3.12), our approach is to mimic the stability 



analysis for the Helmholtz problem (1.4)-(1.6) given in [23l [2ll |35]. It turns out that 
this approach indeed works for scheme ( |3.12[ ) although the analysis is more delicate and 
complicate than that for the differential problem. The key ingredients of our analysis are 
to use a special test function Vh = a - (defined element-wise) with a{x) := x — xq-^ in 

and below) on each element. 



(3.12) and to use the Rellich identity (cf. 

Our first lemma of this section establishes three integral identities which play an 
important role in our analysis. 



Lemma 4.1. Let a(x) 



X 



& E, K ETh and e G ■ Then there hold 



(4.1) 
(4.2) 
(4.3) 



2 / I 1 2 

d\\v\\^2,K) + '2.Re{v,a ■Vv)k = I a ■ nK\v\ , 

JdK 

{d - 2) \\Vv\\%,^. + 2Re {Vv, V(a ■ Vv))^ = / 

JdK 



dv 



'dK 

[a ■ 

[a-Vv] ) - {a -rieiVv} ,[Vv]) 

d-l 



E 

i=i 



a ■ Ti 



dv 
drip. 



a ■ Up 



j dv \ ^ d[ 



dri 



where x^^ denotes the 'point in the star-shaped domain definition for ili (see (2.2) ). Also 



note that in (4.1) and (4.2), we omit the sign ds in the integrals. We shall adopt this 



omission consistently throughout this paper to save space. 

Proof. It is easy to verify by direct computations the following differential identities 
on K: 



div(a(x)) = d, 

diy{avv) = dvv + (a ■ Vv)v + v{a ■ VtJ), 
div(Q!(Vt; ■ Vv)) = {d- 2)Vv ■ Vv + V(a • Vv) -Vv + Vv- V{a ■ Vv). 



(4.1) and (4.2) then follows immediately from integrating the above second and third 



identities over K. 



To prove identity (4.3), from the representations 



Vv 



dv 
dn, 



-n, 



i dv 
'dA 



_ dv 
a ■ Vf = — — a ■ rip 
drip 



d-l 



dv 
d7i' 



ra ■ Ti 
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we have 



dv 
dUf 



|,[a-Vi;]^ -{a-n,{Vv},[Vv])^ 



a ■ Hp 



dv 
drip 



a ■ Up 



a ■ Ti 



dv 
drie 

dv 



d-i 



+ E 



a ■ Ti 



i=i 



dn. 

dv 
dup 



dv 
dup 



a ■ Hp 



dv 
due 

d-l 



dv 

d?e 



( dv 



dv 

d?e 



dv 

d?e 



d[v] 
~dd 



which completes the proof of the lemma. □ 

Remark 4.1. The identity (4.2) can be viewed as a local version of the Rellich identity 
for the Laplacian A (cf IM, Since C E, hence, (|4l|-(|43| also hold for any 

function v = Vh & . 



(4.4) 



Now, taking Vh = Uh in (3.12) yields 

2 



ah{.Uh, Uh) - k \\uh\\ ^2(f^) + ik \\uh\\ i2(r^) = (/, Uh) + {g, Uh)^^ 



Therefore, taking real part and imaginary part of the above equation and using (3.9) and 
(3.10) we get the following lemma. 



Lemma 4.2. LetUhEV^ solve ( |3l2| ). Then 
(4.5) I"/*!?,/, -2Re ^ ^jl^l , K]^ -k'^\\uh\\l2^n)<\if^'>^h) + {g,Uh)rJ, 



(4.6) 



d-l 



70,e II r 111 2 , A.e 



h 



hp 



duh 
drl 



duh 

dup 



LHe) 



+ k \\uh\\\2iY^) < \ {f,Uh) + (fi',Mh)r«| 



From (4.5) and (4.6) we can easily bound \uh\^ ^ and the jumps in terms of \\uh 



In order to get the desired a priori estimates, we now need to derive a reverse inequality 
whose coefficient on the right-hand side can be controlled. Such a reverse inequality, 
which is often difficult to get under practical mesh constraints, and stability estimates for 



scheme (3.12) will be derived next. 

Theorem 4.3. Let Uh € V'^ solve (3.12) and suppose 'Jo e,li e, Pi e > 0. Define 



M{f,g) 



(4.7) 



L2(C) 



+ ||5'||L2(r^). Then there exists a positive constant Cgta such that 



Ml^in) + I hh\\l,h + l{Yl ll^"'»llL2(e) 
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and 



(4.8) Csta < 7 + 77 + 77 max ^ 7— + \ / ^5 ^ ^ 

Proof. Since the proof is long, we divide it into three steps. 

Step 1: A representation identity for \\uh\\]^2(^Qy It follows from (4.1) with v = Uh that 



711 ||2 _ 

Summing over all K E Th yields 



/ a ■nK\uh\'^ -2Re{uh,a ■Vuh)K- 

JdK 



II^/iIIl2(q) = ^ / a ■nK\uh\'^ -2 Re{uh,a ■Vuh)K, 

JdK 



hence, 



(4.9) 2A;^ ||u/i||^2(!^) = /c^ ^ / a ■ nK\uh!\' - {d - 2)k'^ \\uh\\\2^^) - 2k'^Re{uh,Vh), 

Ken -'^^ 

where Vh & E is defined hj Vhlx = ■ "^uhIk for every K E Th. It is easy to check that 
Vh\K is a linear polynomial on K, hence, Vh G V^. Using this Vh as a test function in 
(3.12) and taking the real part of the resulted equation we get 

(4.10) - k'^Re{uh,Vh) = Re{{f,Vh) + {g^Vh)^^ - ah{uh,Vh) - ik {uh,Vh)rJ- 
Now, it follows from (Q, (Q and ( |4.10D that 

2k^ \\uh\\l2(^Q) = ^ « ■ nK\uh\'^ + {d-2) Re[{f,Uh) + {g^Uh)^^^ - a'h{uh,Uh)) 

+ 2 Re((/, Vh) + {g, Vh)^^ - tthiuh, Vh) - ik {uh, Vh)rJ 
= k^ ^ / a ■nK\uh\'^ + {d-2)Re{{f,Uh) + {g,Uh)Y^J 

i^^T. J dK 



(4.11) + 2 Re((/, Vh) + {g, Vh)^^) + 2k Im {uh, Vh)r„ 

- {id-2)\\Vuh\\ 



[vh] 



+ Re(^[uh] , 1^1^ ^ + 2Im {Jo{uh,Vh) + Ji{uh,Vh) + Li{uh,Vh)). 

I 1 2 I 1 2 — 
Using the identity |a| — \b\ = Re (a + 6) (a — 6) we have 

(4.12) ^ / a-nK\uh\^ = 2YRe{a-ne{uh} ,[uh])^ + {a-nn,\uh\^)g^. 

Ken -^^^ e&sl 
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From the Rellich identity (4.2) and noting that (a ■ Ue {Vuh} , [VM/i])e = ■ ne, iVw/il )^ 
for e G Sj^ we get 

=2^Re(a-ne{Vu4,[Vn/,]),+ 



[a ■ Tie, |Vn/,p)^ 



=2 ^ Re (a ■ {VmJ , [Vm^]), + " I'^^'^Oe " 5Z ' I'^^'^He • 



On noting that Uh is piecewise hnear and Vhlx = — xnj ■ Vuhlx, then Vvhlx = Vu^lx- 
Hence, 

(4.14) Im {Ji{uh, Vh) + Li{uh, Vh)) = Ini {Ji{uh, Uh) + Li{uh, Uh)) = 0. 

Plugging ( |4.12D -( |414l ) into ( |4.11[ ) then gives the following representation for H^^h 11^2 (5-2): 

l|Mh||i2(n) = (c^-2)Re((/,Mft) + {g,Uh)TpJ +2Re{{f,Vh) + {g,Vh)rJ 
+ 2k'^ ^ Re (a ■ {uh} , K])^ + k'^ (a ■ uq, luh]"^). 



'on 



(4.15) 



+ 2 J] Re (-{a -He {Vuh} , [Vm/,]), + / | 
+ 25^ ((ci-l)Re/||^|,K]\ +Re/K] 
-2 ^ Re 



[m/i] ) +2lmJo{uh,Vh) 



Step 2: Derivation of a reverse inequality. We bound each terms on the right hand 
side of (4.15). For an edge/face e G S^, let Kf, and K'^ denote the two elements in 
that share e. For an edge/face e G let denote the element in that has e as an 
edge/face and K'^ = ^. We have 



(4.16) 2k^ ^ Re (a ■ {uh} , H), < Ck^ J] /^e"^ h/.|lL^(K.uK^) II Mil 



L2(e) 



< Y Il«^ll22{f^) + ^ E :^|f IIHIli^{e) • 



It is clear that 
(4.17) k'^ {a ■ nn, \uh\'^)g^ = A;^ (a ■ uq, \uh\^)^^ + 



a ■ Up 



< Ck"^ \\uh\\l2(^p^^ + ^ A;^ (a ■ n^, Iw/iP)^ 
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It follows from the star-shaped assumption on Qi that 
(4.18) 



<CkJ2 hh\\L2(e) \\^Uh\\L2(e) ' II VU/, || 



<CA:2|k/.|lL^(r,)-^El|V^^ll'^(e)- 



By (4.3) we obtain 



(4.19) 



duh 

drip 



d-l 



a ■ Ti 



duh 
drip 



d-l 



a ■ Up 



eeSi^ i=l 



K=Ke,K' 



dUh 

dri 



egPlB ^^'^ i=l ^ 



5ri 



Noting that - ^ 

(4.20) 2 ^ (^(rf-i)Re^|^|,K]^ +Re<^K], 

^ E E l|V«.IL.(A')IIK]llL^(e) 



LHe) 



dVh 

dup 



<^K.+^E-^iiMii^ 



70,e 



From (3.3), the inverse inequality and ( |4.6[ ) we get 



(4.21) 2Im (Jo(w/.,^/.)) = 2Im ^ ^ (K] , K])^ 



2Im ^ 



7o^ 
hp 



a ■ Up 



duy, 
dup 



d-l 

+ E 



9ri 



< 



2Im ^ 



70,e / 

— ( a ■ UeUh, ^ — 
Alp \ OUp 



L2(e) 
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duh\ \ d[uh\ 



dri I / dri 



L2(e) 



L2(e) 



hp 



LHe) 



d-l 



L2(e) 



14 



< 



2Im ^ 



70,e 

hp 
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dUh 



drip 



eeS, 



V 7l,e K I /ie 



L2(e) 



/ 70,e / 70,e 



2 , A.e 



duh 
drie 

duh 
dri 



2 



Since D is star-shaped, we have 



(4.22) (" ■ l^'^l'>e + (" ■ l^^'^Oe + ^ ( 



on J 



< _ ^ («-«i?,^V/.r + |VM;,|2-2^|M;,||Vn/,| 



he 



7o.e 7o, 

he he 



7o,e 7o,, 



/ip hp 



2 

L2(e) ' 



Putting (4.15)-(4.22) together we have 



2k' 



L2(n) < (rf- 2) Re + {g,Uh)^^ +2Re{{f,Vh) + (c/,^^h)r«) + y l|Mft|li2(f^) 



70,e /ie 



d-1 



e6£-£0 V"^'^ j=l 



2 



|2 

lL2(e) 



+ 



L2(e) 



1 7o,. 

70,e /ie 



2 

L2(e) 



+ ^E?^?^IK1I' 



/le he 



L2(e) 



/To^ 1 / To, 

ll,ehe \ he 



L\e) 



+ ll,ehe 



duh 
One 
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Ij0,e 



70,e 

hp 



d-l 



iiMiiW) + E 



hp 



duh 



Therefore from ( |4.6[ ), 



< 



1 , 7o,, 



+ C ( + 1 + max + ^ + \lir^ + o 

edS^ V70,e /ie V Pl,e 



7o,. 



1 



+ max 



k"^ + 1 ^ 1 /t^ ^ /t^ ^ 1 



S'tej* 3: Finishing up. It follows from ( |4.5[ ), ( |4.6[ ), (4.8), and the above inequality that 

4A;2 



< 



1| |2 , Cf^i 



where M(/, g) : = 
M{f,gY (cf. (4.6[)) to derive the last inequality. Hence, 



L2(n) 



||5(||i2(r^) and we have used k"^ ||Mft|li2(r^) < k^ 



+ 



\Uh 



which together with (4.6) implies (4.7). The proof is completed. □ 



Remark 4.2. // the penalty parameters are taken as 7i^e = 0, Pi^e = and 7o,e > 0, 



the terms 



dui. 



LHe) 



drl 



L2(e) 



in (4.19) and (4.21) ought he estimated differently 



by using the inverse inequality. This then leads to the following weaker stability estimate: 



(4.23) 



fl , 1 , P + l , l + 7g,e(l + M+7oA , ,,,, 
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Since scheme (3.12) is a linear complex-valued system, an immediate consequence of 



0,/ie > 



the stability estimates is the following well-posedness theorem for (3.12). 

Theorem 4.4. The IPDG method ( |3.12[ ) has a unique solution for k > 
0, 7o,e > 0, a = 1, 7i,e > and (3i^e > 0. 

Remark 4.3. (a) IPDG method (3.12) is well-posed for all wave number k > 
provided that the penalty parameter 7o,e > 0. As a comparison, we recall that [31^ the 
finite element method is well-posed only if mesh size h satisfies a constraint h = 0{k~'^) 
for some r > 1, hence, the existence is only guaranteed for very small mesh size h when 
wave number k is large. 

(b). It is well known that ]3 ^ [7I|, [2Zl |23 123^ symmetric IPDG methods for coercive 
elliptic and parabolic PDEs often require the penalty parameter 7o,e is sufficiently large to 
guarantee the well-posedness of numerical solutions, and the low bound for 7o,e is hard to 
determine and is also problem-dependent. However, this is no issue for scheme (3.12), 



which solves the (indefinite) Helmholtz equation, because it is well-posed for all 7o,e > 0. 



We have the following consequence of Theorem 4^ for quasi- uniform meshes. 
Theorem 4.5. Let h = max/ig. Suppose the mesh is quasi-uniform, that is 



h. 



Assume that 7i^e — 7i > 0, and that 7o,e 
and that 7o,e — 7(f > and [3- 
Ifk>l, then 



> 



{k^hfl^'^l^^ and (3 



> 



i,e ^ lo for e e 4 



D 



i,er:.{h/kfl^^[" force Si, 
where 71 and 7(f are independent of e. 



< 



(4.24) \\uh\\L2in) + ^hhWi^h 
If, furthermore, 7,^ ~ 1 and 71 < k^h, then 



{k^hy/^-f[ 



1/3 



M{f,g). 



(4.25) 



Wh\\L2{n) + 




Mif,g). 



Remark 4.4. It is clear that if k'^h > 1 and 71 and 7^ are chosen properly, say 
7i ^ fci?, 7^* - v^; then 

hhllL^iu) + \ W\\i,h ^ TM{f,g). 



Note that the above estimate is of the same order as the PDF stability estimate given in 
Theorem 2.2. But a large 71 or a small 7,f may cause a large error (cf. Theorem 5.5). 



5. Error analysis. In this section, we derive the error estimates for the solution of 



scheme (3.12). This will be done in two steps. First, we introduce an elliptic projection 
of the PDE solution u and derive error estimates for the projection. We note that such a 
result also has an independent interest. Second, we bound the error between the projection 
and the IPDG solution by making use of the stability results obtained in Section |4j In 
this section, we assume that the mesh is quasi-uniform and 7i^e — 71 > 0. Also, we 
define 



7o := min 7o,e (> 0). 
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5.1. Elliptic projection and its error estimates. For any w & Ed iJp^(fi) fl Hi^^{Q), we 
define its elliptic projection Wh G by 



(5.1) 



ah{wh, Vh) + ik {wh, Vh)r„ = a-h{w, Vh) + ik {w, Vh)i 



In other words, Wh is an IPDG approximation to the solution w of the following (complex- 
valued) Poisson problem: 



dw 



dric 



-Aw 


= F 


in Q, 


w 


= 


on r^i 


+ ikw 




on Tr 



for some given functions F and ip which are determined by w. 

Before estimating the projection error, we state the following continuity and coercivity 
properties for the sesquilinear form ah{-, ■)• Since they follow easily from (3.1)-(3.10), so 
we omit their proofs to save space. 

Lemma 5.1. For any v E E and w E E r\ Hy^{VL), the mesh- dependent sesquilinear 
form ah{-,-) satisfies 



(5.2) 



\ah{v,w)\, \ah{w,v)\ < \\v\\^^,^ \\\w\\\i,h- 



In addition, for any < e < 1, there exists a positive constant such that 
(5.3) Re ah{vh, Vh) + [1 - e + — ) Im ah{vh, Vh) > (1 - £ 



as above. Then (5.1) immediately implies the following Galerkin orthogonality: 

(5.4) ah{u-Uh,Vh) + \k{u-Uh,Vh)^^ = ^ ^vneV^. 

Lemma 5.2. Suppose problem (1.4)-(1.5) is H'^ -regular. Then there hold the following 
estimates: 

(5.5) 
(5.6) 



I" ^ ^hWi^h + i^^f^ 11^ - ^h\\L2(TR) ~ ^(^ + 7i + kh)^kh, 
\\u - M/i||i2(Q) < A(A + 7i + kh)kh^, 



70 



where A := 1 

Proof. Let Hh be the Pi-conforming finite element interpolation of u on the mesh Th. 
Then Uh & E D H^^^Q) and satisfies the following estimates (cf. [TH [T8]): 



(5.7) 

(5.8) 

(5.9) 

Let r]h 
(5.10) 



1^ '"^llL2{n) ~ ^^I'"l/f2(f7)5 



I" - ^h\\\i,h ~ 1 +71 + 



7o 



h\u\ 



(A + 7i)"/i| 



u 



\u-Uh\\L2(^rR) 



Uh — Uh- From rjh + u — Uh = u — Uh and (5.4) 



a-hiVh, Vh) + ik {r]h, 'nh)^ = (^h{u - Uh, Vh) + ik {u - Uh, Vh)ri, ■ 
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Take £ = | in (5.3) and assume without loss of generality that ci > ^. It follows from 
Q and dsloj ) that 

1 2 /I 

2 ' \2 70/ 

/I Cl X 

= Re(a,,(u - Uh,Vh) + ik {u - Uh.rjh)^^ ^ V2 —)f^{Vh,Vh)rji 



+ 



1 Cl N 

- + ^ j Im [ah{u - Uh, r]h) + ik {u - Uh, Vh)rji) 



<CA^||?7/,||i^;^ |||n - UhWli^h + k\\u- Uh\ 



A/c 2 

4 WVhWL^Fn)- 



Therefore, it follows from (5.8), (5.9) and (2.2) that 

|2 



(5.11: 



VhWi^h + \\Vh\\L2(rR) ~^ 111^ - ^hlWi^h + h - ^h\v^^^Vn) 

<X^eh^{X + ^i + kh). 



which together with the fact that u — Uh = u — iih ~ yields (5.5). 

To show (5.6), we use the Nitsche's duality argument (cf. [m HE])- Consider the 
following auxiliary problem: 



(5.12) 

dw 

It can be shown that w satisfies 
(5.13) \w\ 



-Aw = u — Uh 
w = 

ikw = 



in Q, 

on r^, 

on Tr. 



< 



Let Wh be the Pi-conforming finite element interpolation of w on Th. Testing the conju- 
gated (5.12) hy u — Uh and using (5.4) we get 

11^ - ^/i|lL2(n) = -{u-Uh,Aw) = ah{u-Uh,w) + ik {u-Uh,w)^^ 
= ah{u -Uh,w - Wh) + ik{u- Uh, w - w/i)rp 



< \\U - Uh\\i 



\w — 



< 



\u 



^hWi^h (-^ + li)^h\w 



+ k\\u- Uh\ 



h'2\w\ 



which together with (5.5) and (5.13) gives ( |5.6[ ). The proof is completed. □ 

5.2. Error estimates for scheme ( |3.12[ ). In this subsection we shall derive error estimates 
for scheme (3.12). This will be done by exploiting the linearity of the Helmholtz equation 
and making use of the stability estimates derived in Theorem |4.3 and the projection error 
estimates established in Lemma 15.21 

Let u and Uh denote the solution of (1.4)-(1.5) and that of (3.12), respectively. Assume 
that u G H'^{Q). Then (3.11) holds for v = Vh & V^. Define the error function eh '■= u — Uh- 
Subtracting (3.12) from (3.11) yields the following error equation: 

yvh e V\ 



(5.14) 



ah{eh, Vh) - k {eh, Vh) + \k{eh, Vh) 
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Let Uh be the elliptic projection of u as defined in the previous subsection. Write et = rj—^ 
with 



V 



u 



i := Uh - Uh. 



From (5.14) and (5.4) we get 
(5.15) 



ah{^, Vh) - k Vh) + i/c(^, Vh)TR = ahiv, ^ft) - ^ iv, ^h) + Uh)rR 

= ~k\ri,Vh) \/vheV\ 



The above equation implies that ^ G is the solution of scheme (3.12) with source terms 
/ = —k'^T] and g = Then an application of Theorem 4^ and Lemma 5^ immediately 
gives the following lemma. 

Lemma 5.3. ^ = Uh — Uh satisfies the following estimate: 



(5.16) 



llell 



L2(n) 



lieiU<ataA(A + 7i + A;/i)fc3/i2 



70 



where Cste, is defined in (4.8) and A = 1 + 

We are ready to state our error estimate results for scheme (3.12), which follows from 
Lemma |5.3[ Lemma |5.2| and an application of the triangle inequality. 



Theorem 5.4. Let u and Uh denote the solutions of (1.4)-(1.5) and (3.12), respec- 
tively. Assume that u G H'^iVL). Then there exist two positive constants Ci and C2 such 
that the following error estimates hold. 



(5.17) 
(5.18) 



\u 



\u 



- w^lli,^ < A(A + 71 + kh) [CM + C^C^.^k^h''), 
Mh|li2(f^) < A(A + 71 + kh) {Cikh^ + C2C,te.k^h^), 

1 



where Cgta. is defined in TheoremU.o] and A = 1 + 



70 



From Theorem 4.5 and the definition of Csta (cf. (4.8)) we obtain the following esti- 
mates. 

Theorem 5.5. Assume that < 71 < k'^h, and that 7o,e — {k'^hy^^'yl^^ and (3i^e > 
{h/ky^^'yl^^ for e G S^, and that 7o,e — 1 and /5i_e > 1 for e & Sj^ . If k > 1, then there 
exist two positive constants Ci and C2 such that the following error estimates hold. 



(5.19) 
(5.20) 



\\u-Uh\\^^f^ < A(A + 7i + {Cikh + C2(^ + 
\u-Uh\\L2(a) < \{\ + -ii + kh){Cikh'^ + C2(^ ' 



{k^hy/^-il''^ 
1 



k^h'' 



{k^hy/^-f[ 



1/3 



Remark 5.1. (a) The estimates in (5.17)-(5.20) are so-called preasymptotic error 
estimates (i.e. for the mesh in the regime k^h > 1). In fact, the estimates hold for any 



h > 0. We recall that [37^ the preasymptotic error estimates for the finite element method 
solution was only proved in the 1-d case provided that kh < 1. 

(b) The second term on the right-hand side of the first inequality is pollution term for 

h-Uh\\i^h- 

(c) If kh <1, then under the assumption of Theorem 5.5 we have 
(5.21) \\u - UhW^^h < Cikh + C2k^'^h^'^ 



20 



X. FENG AND H. WU 



for some constants Ci and C2 which depend on 71. Numerical tests in the next section 
suggest that \\u — Uh\\if^ may have a better bound Cikh + C2k^h'^ and it is possible to tune 
the penalty parameters to significantly reduce the pollution error. We note that in the case 
k'^h is sufficiently small, optimal order (with respect to h) error estimate in the broken 
H^-norm can be derived by using the Schatz argument as done in 13^ [13 \27^ . 
(d) Inequality (5.16) shows that \\uh ~'^h\\ih ^iT'joys a superconvergence. 

6. Numerical experiments. Throughout this section, we consider the following two- 
dimensional Hclmholtz problem: 

sin(A;r) 



(6.1) 
(6.2) 



-Am 

du 

dnn. 



-k'u = f :-- 
+ iku = g 



in 



on r 



R 



dQ. 



Here f2 is the unit regular hexagon with center (0,0) (cf. Figure 6.1) and g is so chosen 
that the exact solution is 



(6.3) 



u 



cos{kr) 



cos k + i sin k 



-Mkr) 



k k{Jo{k) + Uiik)) 
in polar coordinates, where Ju{z) are Bessel functions of the first kind. 




Fig. 6.1. Geometry (left) and a sample mesh Ti/^ that consists of congruent and equilateral triangles of size h = 1/7 
(right) for Example 1. 

For any positive integer m, let Ti/m denote the regular triangulation that consists of 
6m^ congruent and equilateral triangles of size h = 1/m. See Figure 6.1 (right) for a 
sample triangulation Ti/j. 

6.1. Stability. In this subsection, we use the following penalty parameters for the 
symmetric IPDG method (cf. (3.12)) according to Theorem 4.5| (or 5.5): 

0.1, 7o,e = {k'hr/'iy\ 



(6.4) 



71. 



and (3 



l,e 



1. 



Given a triangulation T^, let m™'^ be the Pi-conforming finite element approximation 
of the problem (6.1)-(6.2). Recall that Uh denotes the IPDG solution. Figure 6.2 plots 



DG METHODS FOR THE HELMHOLTZ EQUATION 



21 



the if^-seminorm of the IPDG solution the iJ^-seminorm of the finite element 

for h = 0.05 and 0.005, respectively, and the if^-seminorm of the 
for k = 1, - ■ ■ , 230. It is shown that 



solution 

exact solution \u 



(6.5) 



\u 



\u 



FEM| 
h 



< 1. 



We notice that the stability estimate \\uh\ 
stability estimates are better than those given by Theorem |4.5 



< 1 implies that 



ML2(n) ^ lA- These 




250 



Fig. 6.2. Ilufelli (solid), ["^^^l^fUfj) (dashed) for h = 0.05 and 0.005, respectively. The dotted line gives the 
H^-seminorm of the exact solution \^\uii^Qy 



6.2. Error of the finite element interpolation. Given a triangulation Th, let Uh be the 



Pi-conforming finite element interpolation of u on 7^. Consider in Figure [673] log- log plots 
of the relative error e{h,k) := |m — u/il^j^/lul^ of the finite element interpolation in H^- 
seminorm for different k versus 1/h. Similar to the 1-D case [37], AH error curves decay 
with constant slope of —1. Note that the error stays at around 100% on coarse mesh and 
starts to decrease at a certain mesh size. We are interested in the mesh size where the 
descent starts. Similar to "the critical number of degrees of freedom" introduced in |37j . 
we introduce the following definition of critical mesh size. 

Definition 6.1. Define — for any fixed k and f — the critical mesh size as maximum 
mesh size H{k, /) for which 

1. e{h, k) <1 for h < H{k, f), and 

2. e{h,k:) -^0 as 

Recall that the critical mesh size for the one dimensional case is one half of the 
wavelength, that is, | Since the solution u is axial symmetric (cf. (6.3)) and the 

trace along any direction may be resolved by a mesh with mesh size less than |, the 
critical mesh size for the finite element interpolation should be greater than or equal to 



Figure [674] plots the reciprocal of the critical mesh size for the finite element interpolation 
computed for all integer k from 1 to 230 and the line passing through the origin has slope 



— ^7=. It shows that 



(6.6) 



the critical mesh size for Uh 



73 
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Fig. 6.3. Relative error of the finite element interpolation in -seminorm for k = 5, k = 10, k = 50, and k = 100. 
The dotted line gives the reference line with slope —1. 




250 



Fig. 6.4. Reciprocal of the critical mesh size of the relative error of the finite element interpolation in -seminorm 



computed for k = 1, 



, 230. The dotted line gives the line through the origin with slope 



Figure 6^ also shows that the error of the finite element interpolation is controlled 
by the magnitude kh. For illustration, the points that are computed from kh = 0.25 
are connected. The connecting line does neither increase nor decrease significantly with 
the change of k. For more detailed observation, the relative errors of the finite element 
interpolations, computed for all integer k from 1 to 230 for kh = 1 and kh = 0.5, are 

1 stays around 0.247 and the error for kh = 0.25 
0.5 which verifies that the relative error of 
0{kh). 



plotted in Figure |6.5[ The error for kh 
stays around 0.124. Note that 0.124/0.247 
the finite element interpolation in if^-seminorm satisfies e{h, k) 



6.3. Error of the DG solution. From Theorem 5.4, the stability estimates in Subsec- 



tion 

by 

(6.7) 



6.1 suggest that the error of the IPDG solution in if^-seminorm could be bounded 



\u~Uh\^^< {l + kh){Cikh + C2k^h^) 
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Fig. 6.5. Relative errors of the finite element interpolations in H^- 
size h determined by kh = 1 and kh = 0.5, respectively. 



leminorm computed for k = 1, 



, 230 with mesh 



for some constants Ci and C2. The second term on the right hand side is the so-called 
pollution error. We now present numerical results which verify the above error bound. 



In Figure 6.6, the relative error of the IPDG solution with parameters given by (6.4) 
and the relative error of the finite element interpolation are displayed in one plot. The 
relative error of the IPDG solution stays around 100% before a critical mesh size is reached, 
then decays slowly on a range increasing with k, and then decays at a rate greater than 

— 1 in the log-log scale but converges as fast as the finite element interpolation (with slope 

— 1) for small h. The relative error grows with k along line kh = 0.25. Unlike the error of 
the finite element interpolation, the error of the IPDG solution is not controlled by the 
magnitude of kh — see also Figure |6]7 



10 



10 



o 
Si 

g 10" 



10 



10 



10 



10 



iiti 



Fig. 6.6. Relative error of the IPDG solution with parameters given by | |6.4| (solid) and relative error of the finite 
element interpolation (dotted) in -seminorm for k = 5, k = 10, k = 50, and k = 100, respectively. 



Figure 6.8 plots the relative error of the IPDG solution with parameters given by (6.4) 
for k = 52/^0^/3, . . . , 500^/3 and h determined by k^h"^ = 1. The error does not increase 
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250 



Fig. 6.7. Relative errors of the IPDG solutions with parameters given by | |6.4| l in -seminorm computed for k 
1, • • ■ ,230 with mesh size h determined by kh = 1 and kh = 0.5, respectively. 



with respect to k which verifies (6.7). 




Fig. 6.8. Relative errors of the IPDG solutions with parameters given by \&A\ in -seminorm computed for k 
5^/^, 10^/^, ■ ■ ■ , 500^/^ with mesh size h determined by k^h"^ = 1. 



Figure |6.9| plots the reciprocal of the critical mesh size for the IPDG solution with 

il 

and -, respectively. It shows that 



parameters given by (|6.4|) computed for all integer k from 1 to 230 and the lines passing 

1 



through the origin have slopes ^ 
(6.8) 



the critical mesh size for Uh 



l.SSvr 



with two exceptions but they are still less than |. It is interesting that the dependence on 
is essentially linear. We consider the IPDG solution with parameters given by (6.4) 
for k = 100 on the mesh with mesh size h = 1/60. The relative error in if^-seminorm 



is about 0.9898. Figure 6.10 presents the surface plots of the interpolation (left) and the 
IPDG solution (right). It is shown that the IPDG solution has a correct shape although 
its amplitude is not very accurate. 
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250 



Fig. 6.9. Reciprocal of the critical mesh size of the relative error of the IPDG solution with parameters given by ( |6.4[ l 
-seminorm computed for k = 1, - ■ ■ , 230. The dotted lines give the lines through the origin with slopes y^s^ '^"'^ ^ • 
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Fig. 6.10. Surface plots of the interpolation (left) and the IPDG solution (right) with parameters given by l |6.4[ l for 
100 on the mesh with mesh size h = 1/60. 



6.4. Sensitivity of the error bounds with respect to penalty parameters. In this subsec- 
tion, we examine the sensitivity of the error of the IPDG solution in if ^-seminorm with 
respect to the parameters 7o,e, Pi,e, and 71^6, respectively. 

First, we examine the sensitivity of the error in 7o,e- To the end, for k = 5 and 
k = 50, respectively, we fix ■ji^^ = 0.1 and Pi^e = 1, and compute the IPDG solution with 
the following 70,, 



7o,( 



{k'hyHH (see dp)), 7o„ 



70,e 



0.01, and 7o,e = 100. 



Figure [6.11| plots the relative error of the IPDG solution for each run. We observe that 
the error in the if ^-seminorm is not sensitive with respect to the parameter 7o,e- It is 
clear that 7o,e affects the continuity of the solution. The larger the parameter 7o,e, the 
more continuous the IPDG solution. 



Secondly, we test the sensitivity of the error in /?i e. Figure 6.12 plots the relative 
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Fig. 6.11. Relative error of the IPDG solution with parameters 71 e = 0.1,/3i,e = 1; o.'<^d, each of following 70, e- 

70, e = (fc^/i)'^'''^7i{f (dotted), 70, e = 1 (solid), 70, e = 0.01 (dashed), and 70, e = 100 (dashdot) in the -seminorm for 
k = 5 and k = 50, respectively. 

error in the if^-seminorm of the IPDG solution with parameters 7o,e = l;7i,e = 0.1, and 
each of the following Pi^g-. /3i,e = 0, 1, 100, for = 5 and k = 50, respectively. Again, we 
observe that the error in the i7^-seminorm is not sensitive with respect to the parameter 

Pl,e. 
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Fig. 6.12. Relative error in the -seminorm of the IPDG solution with parameters 70, e = li7l,e = 0.1, and each of 
the following Pi,e: /3i,e = (dotted), /3i,e = 1 (solid), and /3i,e = 100 (dashed) for k = 5 and k = 50, respectively. 



Finally, we examine the sensitivity of the error in 71 g. To the end, we fix 7o_e = 1 and 
Pi^e = 1 and compute the IPDG soluti on wit h the following ji^e- 7i,e = 0,0.01,0.1, 1 for 
k = 5 and k = 50, respectively. Figure 6.13 plots the relative error in the if^-seminorm 
of the IPDG solution for each run. We observe that the error has a similar behavior as 
the error of the finite element solution for small value 71^, 
the solution is more stable for larger 71 g, 
large (absolute) error. 



0,0.01 (cf. Figure [6l6| ) and 
but a large 71 _e, say 71 = 1, may result in a 
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Fig. 6.13. Relative error in the H^-seminorm of the IPDG solution with parameters 70, e = li/9l,e = 1, and each of 
the following 71, e-' 71, e = 0.1 (solid), -yi^e = (dotted), 71, e = 0.01 (dashed), and ■yi^e = 1 (dashdot) for k = 5 and k = 50, 
respectively. 



6.5. Reduction of the pollution effect. One advantage of the IPDG method is that it 
contains several parameters which can be tuned for a particular purpose. In [2], it is 
shown that it is possible to reduce the pollution error of the IPDG method by choosing 
appropriate parameters a and i7o,e- Recall that all choice of a but one lead to non- 
symmetric formulations. In this subsection, we shall show that appropriate choice of the 
parameter 7i_e can significantly reduce the pollution error of the symmetric IPDG method 
(with (7 = 1). We use the following parameters: 



(6.9) 



171,- 



-0.07 + O.Oli, 



7o,, 



100, 



and j3 



l,e 



1. 



We remark that i'^i,e is simply chosen from the set {0.01(p + gi), — 50 <p,q < 50} to 
minimize the relative error of the IPDG solution in if^-seminorm with 7o,e = 100 and 
Pi^e = 1 for wave number k = 50 and mesh size h = 1/20. 



In Figure 6.14, the relative error of the IPDG solution with parameters given by (6.9) 



and the relative error of the finite element interpolation are displayed in one plot. The 



IPDG method with parameters given by (6.9 ) is much better than the IPDG method using 
parameters given by (6.4) (cf. Figure 6.6). The relative error does neither increase nor 
decrease significantly with the change of k along line kh = 0.25 for k < 100. But this does 
not mean that the pollution error has been eliminated. For more detailed observation. 



the relative errors of the IPDG solution with parameters given by (6.9), computed for all 
integer k from 1 to 230 for kh = 1 and kh = 0.5, are plotted in Figure |6.15[ I t is shown 
that the pollution error is reduced significantly (cf. Figure 6.7 and Figure 6.5). 



6.6. Comparison between the IPDG solution and the finite element solution. We have 
shown the flexibility and performance of the IPDG method in previous subsections. In 
this subsection, we give a comparison between the IPDG method and the finite element 
method. One disadvantage of the IPDG method compared to the finite element method 
is that the linear system of the IPDG discretization involves more number of degrees of 
freedom than that of finite element discretization on the same mesh. In two dimensional 
case it is about three times more. So in the asymptotic range, the IPDG method is less 
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Fig. 6.14. Relative error of the IPDG solution with parameters given by 
element interpolation (dotted) in -seminorm for k = 5,k = 10, k = 50, and k 



9|l (solid) and relative error of the finite 
100, respectively. 




250 



Fig. 6.15. Relative errors of the IPDG solution with parameters given by \6.9\ in -seminorm computed for k 
1, ■ ■ ■ , 230 with mesh size h determined by kh = 1 and kh = 0.5, respectively. 



effective in terms of number of degrees of freedom. We shall show that, for Problem 



(6.1 )-(6.2), the IPDG solution is more stable than the finite solution for large h, and it is 



possible to choose appropriate parameters such that the IPDG method is more effective 
than the finite element method in preasymptotic range even in terms of number of degrees 
of freedom. 



In Figure |6.16[ the relative error of the finite element solution and the relative error of 
the finite element interpolation are displayed in one plot. The relative error of the finite 
element solution first oscillates around 100%, then decays at a rate greater than —1 in 
the log-log scale but converges as fast as the finite element interpolation (with slope —1) 
for small h. The relative error grows with k along line kh = 0.25. The error of the finite 
element solution is not controlled by the magnitude of kh — see also Figure |6.17 



Figure 6.18 plots the reciprocal of the critical mesh size for the finite element solution 
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Fig. 6.17. Relative errors of the finite element solutions in -seminorm computed for k = 1, - ■ ■ , 230 with mesh size 
h determined by kh = 1 and kh = 0.5, respectively. 



computed for all k from 1 to 230 and the curve m = ^Jk^/AS. It is shown that 



(6.10) 



the critical mesh size for uj™ 




We can see that the IPDG solution is more stable than the finite element solution. 



For more detailed comparison, we consider the problem (6.1)-(6.2) with wave number 
k = 100. The traces of the IPDG solutions with parameters given by (6.9) and the finite 
element solutions in the a;2;-plane for mesh sizes h = 1/50, 1/120, and 1/200, and the trace 



of the exact solution in the x2;-plane, are plotted in Figure 6.19 The shape of the IPDG 
solution is roughly same as that of the exact solution for h = 1/50,. They matches very 
well for h = 1/120 and even better for h = 1/200. While the finite element solution has a 
wrong shape near the origin for h = 1/50 and h = 1/120 and only has a correct shape for 
h = 1/200. The phase error appears in all the three cases for the finite element solution. 
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Fig. 6.18. Reciprocal of the critical mesh size of the relative error of the finite element solution in -seminorm 
computed for k = 1, - ■ ■ , 230. The dotted line gives the curve m = ^ /AS. 



Table 6.1 shows the numbers of degrees of freedom needed for 30% relative errors 
in i7^-seminorm for the finite element interpolation, the IPDG solution with parameters 



given by (6.9), and the finite element solution, respectively. The finite element method 
needs less DOFs when = 10 and /c = 50 than the IPDG method does, but the situation 
reverses when = 100 and k = 200. 



k 


10 


50 


100 


200 


Interpolation 


217 (1/8) 


5,167 (1/41) 


20,419 (1/82) 


81,181 (1/164) 


IPDG 


1,152 (1/8) 


38,088 (1/46) 


217,800 (1/110) 


1,431,432 (1/282) 


FEM 


397 (1/11) 


30,301 (1/100) 


229,357 (1/276) 


1,804,201 (1/775) 



Table 6.1 

Numbers of degrees of freedom needed for 30% relative errors in -seminorm for the finite element interpolation, 
the IPDG solution with parameters given by l |6.9[ l, and the finite element solution, respectively. The fractions in the 
parentheses give the corresponding mesh sizes. 
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